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RADaR  cross  section  evaluation  of  .arbitrary  cylinders 

INTRODUCTION 

The  objective  of  this  program  has  been  to  develop  a  2D  code  and  numeri¬ 
cal  technique  for  determining  the  Radar  Cross  Section  (RCS)  and  other 
electromagnetic  scattering  behavior  of  an  infinite  cylinder  of  arbitrary 
cross-sectional  shape  and  material  composition.  Both  monostatic  and  bis¬ 
tatic  RCS's  need  be  obtained. 

The  technique  we  have  used  to  meet  this  objective  differs  greatly  from 
conventional  methods  of  RCS  evaluation.  Rather  than  working  in  the  fre¬ 
quency  domain,  we  have  illuminated  the  target  by  an  Electromagnetic  Pulse 
(FNP)  containing  significant  energy  over  all  frequencies  of  incerest.  We 
have  then  solved  the  near- field  problem  by  Time -Domain  Finite  Differencing 
(TDFD)  ,,  and  Fourier  transformed  the  result  to  obtain  the  desired  RCS  as  a 
function  of  frequency.  A  near-field  to  far- field  transformation  is  also 
incorporated  into  our  Fourier  procedure. 

Our  technique  has  been  tested  with  targets  up  to  5  m  across  m  obtain 
results  from  50  to  500  MHz.  (These  parameters  correspond  to  a  maximum  of  8 
free-space  wavelengths.)  There  is  nc  reason  we  could  not  obtain  information 
up  to  1  GHz,  although  we  have  not  actually  done  this  to  date  in  the  interest 
of  conserving  computer  resources. 

We  have  run  test  cases  using  generalized  lossy  inhomogeneous  media, 
where  the  real  and  imaginary  parts  of  the  permittivity  or  permeability  are 
frequency  and  position  dependent.  Additionally,  the  code  is  generalized  to 
treat  materials  with  anisotropy  in  the  plane  perpendicular  to  the  cylinder 
axis.  The  code  can  also  treat  generalized  material  discontinuities  includ¬ 
ing  resistive,  capacitive  and  conductive  cards. 

Additionally,  we  have  obtained  canonical  solutions  to  several  problems 
with  known  frequency- domain  solutions  in  order  to  evaluate  the  accuracy  of 
our  code.  These  problems  include  a  circular  cylinder  composed  of  an  ar¬ 
bitrary  number  of  concentric  cylindrical  shells,  where  the  electrical 
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parameters  of  each  shell  are  arbitrary  (but  isotropic) .  We  have  also  com¬ 
puted  the  scattering  by  a  thin  dielectric  strip  of  infinite  extent,  using 
the  technique  of  Richmond.^  Code-code  comparisons  between  the  canonical  and 
TDFD  results  for  these  two  problems  will  be  presented  later. 

Additionally,  we  have  derived  the  equations  for  scattering  by  a  per¬ 
fect,  circular  dielectric  cylinder  with  anisotropy  in  the  plane 
perpendicular  to  the  cylinder  axis.  At  present,  we  have  not  had  time  to 
code  up  this  canonical  solution. 

Also,  we  have  worked  out  the  Mathieu- function  expansion  for  scattering 
off  a  perfect,  elliptic  dielectric  cylinder.  This  solution  has  been  coded, 
and  gives  answers  which,  in  the  perfect-conductor  limit,  agree  with  known 
results.  We  have,  as  of  yet,  not  run  a  comparative  TDFD  case.  This 
Mathieu- function  expansion  has  been  generalized  to  treat  the  case  of  a  two- 
media  elliptical  cylinder  with  the  outer  medium  confocally  coating  the 
inner.  For  our  two -medium  case,  ...he  inner  medium  may  he  perfectly 
conducting.  The  two -medium  solution  hoc  not  yet  been  coded,  nor  r.  compara¬ 
tive  TDFD  run  made. 

In  the  following  sections  of  this  report,  we  will  first  document  the 
TDFD  code,,  Then  we  shall  describe  the  canonical  solutions  and  present  code¬ 
code  comparisons  for  those  cases  where  comparisons  have  been  performed. 
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ARCHITECTURE  OF  THE  TIME-DOMAIN  FINITE-DIFFERENCE  RCS  CODE 


Equations  to  be  Solved 


Ue  assume  in  this  study  that  eiectromagnetically  linear  conditions 
prevail.  Then  the  total  field  can  be  separated  into  an  incident  field 
(which  would  be  the  field  in  the  absence  of  the  scatterer)  and  a  scattered 
field  (the  field  modification  caused  by  the  scatterer' s  presence): 
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+  (£: 


scat 
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Under  the  linearity  assumption,  both  the  incident  and  the  scattered  fields 
individually  satisfy  Maxwell's  equations. 

Tf  some  background  dissipation  (o,  ,  cr,  )  is  present,  the  incident  fields 
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In  the  presence  of  an  anisotropic  scatterer  with  frequency  dependent 
properties,  the  total  fields  conform  to 
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Subtraction  of  the  first  pair  of  equations  from  the  second  leaves  us  with 
the  version  of  Maxwell's  equations  obeyed  by  the  scattered  fields: 
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It  is  convenient  to  represent  the  inhomogeneous  parts  (or  incident 
parts)  of  eqs .  (6)  and  (7)  as 


[JT*]  =  -  [cr*]sHS]  +  [J*]  +  [V  X  ESCat] 


u 
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Thus,  [H  ]  and  [E  ]  become 

[HS]  -  -  [r*][HinC]  -  [R*][HitlC]  -  J  [A*(t-t')][HinC(t')]dt'  (10) 

-00 

[ES]  -  -  [r)[EinC]  -  (R][EinC]  -  J*  [A(t-t')][Einc(t')]dt'  (11) 

-CO 

In  the  last  pair  of  equations,  we  have  defined 

[r*]  -  Mol  (12) 

[R*3  -  [4]'1[4-  <3  (13) 
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[A*]  -  [oJfV*]  (14) 


Method  of  Solution 


For  now,  we  shall  only  concern  ourselves  with  the  2D  TM  case.  Thus, 

E  ,  E  and  H  will  be  present.  Equation  (18)  reduces  to  a  scalar  equation, 
X  y  z 

but  (19)  remains  a  two -component  vector  equation. 

Figure  1  illustrates  a  typical  TDFD  2D  unit  cell,  and  shows  wiere  the 
three  field  components  associated  with  that  cell  are  located. 

In  Appendix  1,  we  will  derive  a  technique  for  solving  this  equation 
system  using  first-order  exponential  differencing.  The  result  of  this 
appendix  is  that,  if  we  omit  frequency  dependence,  HZ(I,J)  is  advanced 
according  to 
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and  [E(I,J)]  is  advanced  according  to 
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(The  complication  of  frequency  dependence  will  be  considered  later.) 


As  we  have  said,  [a0]  is  permitted  to  be  anisotropic.  In  the  actual 
code,  it  Is  represented  by  a  total  of  five  arrays: 
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SGX(I,J)  is  the  bulk  value  characterizing  the  xx  component  of  the 
conductivity  tensor  at  cell  (1,J) 

SGY(I,J)  is  the  bulk  value  characterizing  the  yy  component  of  the 
conductivity  tensor  at  cell  (I,J) 

SGXY(I,J)  is  the  bulk  value  characterizing  the  xy  and  yx  components  of 
the  conductivity  tensor  at  cell  (I.J)  (gyrotropic  materials  are  not 
permitted  in  the  present  code) 

SGCX(I.J)  is  the  surface  conductivity  in  the  x  direction  on  the  y- 
facing  surface  of  cell  (I,J) 

SGCY(I , J)  is  the  surface  conductivity  in  the  y  direction  of  the  x- 
facing  surface  of  cell  (I,J) 

Thus,  the  actual  conductivity  seen  by  Ex  at  cell  (I,J)  is  given  by 


CT°(I*J)xx  "  <SGX(I,J-1)  +  SGX(I,J))/2  +  SGCX(I,J)  (22) 

“  (SGXY(I.J-l)  +  SGXY(I , J) )/2  (23) 

x.y 

"  (SGXY(I.J-l)'1  +  SGXYd.J)'1)*1  *  2  (24; 

yx 

^o(I.J)yy  -  (SGY(I.J-l)'1  +  SGYd.J)'1)'1  •  2  (25) 


This  arrangement 
parallel,  while 


occurs  because  J  sees  the  xx  and  xy  conductivities  in 
x  J 

sees  the  yx  and  yy  conductivities  in  series. 


The  conductivity  matrix  for  E^. 


at  cell  (I,J)  is 


analogously  described 
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°r°(I’J)xx  “  (SGXd-l.J)"1  +  SGX(I , J) _1)  1  •  2  (26) 

<r0(I.J)vv  “  (SGXY(I*1 , J)  1  +  SGXYd.J)*1)’1  •  2  (27) 

XJ 

o0(I,J)  -  (SGXY(I-l.J)  +  SG5CY ( i , J ) ) /2  (28) 

cTod.J)^  -  (SGY(I-l.J)  +  SGY(I , J) )/2  +  SGCY(I,J)  (29) 


Note  that,  although  the  physical  conductivity  tensor  is  symmetric  at 

each  cell  (SGXY(I,J)  —  SGYX(I,J)),  the  mathematical  conductivity  just 

described  is  not  symmetric  (<r0(I,J)  r*  cr0(I,J)  ). 

xy  yx 

The  dielectric  properties  of  the  scatterer  are  represented  by  five 
analogous  arrays,  EPX,  EPY,  EPXY,  EPCX  and  EPCY.  These  are  combined  in  the 
same  way  to  form  the  mathematical  permittivities  e^d.J)^  at  the  Ex  and 
evaluation  point?  of  cell  (I,J). 


Due  to  the  anisotropic  cross-terms,  it  is  necessary  to  know  at  the 
Ex  evaluation  points.  This  is  done  by  simple  linear  interpolation, 


EY(I,J)x  -  (EY(I.J)  +  EY(I+1,J)  +  EY(I , J-l)  +  EY(I+1 , J-l) )/4  (30) 

Ex  at  the  evaluation  points,  EX(I,J)^  is  obtained  the  same  way.  The 

matrix  difference  equation  (21)  is  then  solved  twice  at  each  cell  and  each 
time  step,  once  centered  at  and  to  advance  EX(I,J),  and  once  centered  at  and 
to  advance  EY( I , J ) . 

The  following  notation  is  also  used; 
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QXX(I.J)  is.  the  (1,  1)  component  of  e 
points 


-lO  t^o)At 


evaluated  at  the  E 


QXY(I,J)  is  the  (1,  2)  component  of  e 
points 


-[e  fVolAt 


evaluated  at  the  E 


-U^r1  [<*<)]  At 

QYX(I,J)  is  the  (2,  1)  component  of  e  evaluated  at  the  E^ 

points 

-UJ-^olAt 

QYY(,i,J)  is  the  (2,  2)  component  of  e  evaluated  at  the 

points 


,-l, 


*[ej  [t^0 1  At  . 

SXX(I,J),  SXY(I, J) ,  SYX(I , J)  and  SYY(I,J)  are  ([I]  -  e  )[o0]  L 

correspondingly  located  and  defined. 


Additionally,  if  we  refer  back  co  eqs,  (15)  ■  (17), 


TAUXX(I.J),  TAUXY ( I , J ) ,  TAdYX(I,J)  and  TAUYY(I,J)  are  t^o ) _1 C •  e0] 
analogously  located  and  defined,  and 

RXX(I ,  J) ,  RXY(I ,  J) ,  RYX(I,J>  and  RYY(I,J)  are  [a0]'V0  -  ob] 

analogously  located  and  defined. 


g 

It  is  necessary  to  evaluate  ooth  components  of  [E  j  as  given  in  eq. 

(11)  at  both  E  and  E  points  in  each  cell.  The  above  conventions  indicate 
x  y 

how  to  combine  the  aQ  and  tensors  for  an  inhomogeneous  scatterer  so  this 

complete  set  of  evaluations  may  be  achieved.  In  particular,  we  denote  E^ 

S  c  S 

as  E  evaluated  at,  the  E  mesh  point  and  E  as  E  evaluated  at  the  E  mesh 
y  x  r  yx  x  y 

point. 
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Thus ,  £^(1 , J)  is,  from  eq.  (11)  with  frequency  dependence  not  yet 


included, 


4(I,J)  -  -  TAUXX(I ,  J)E^C(I,  J)  -  TAUXY(IJJ)E^C(I>J) 


-  RXX(I,J)E^C(I,J)  -  RXY(I,J)E*"C(I,J) 
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where  E  (I,J)  is  the  analytically  specified  E  field  evaluated  at  E 
XX  t  X  «  X 

mesh  points,  and  E^C(I,J)  is  he  analytically  specified  Ey110  field 

evaluated  at  E  mesh  points.  Additionally,  ES  (I,J)  is 
x  yy 


1  S 


Eyy(l , J)  -  -  TAUYY(I , J)E^C(I, J)  -  TAIWX(I,J)E^C(I,J) 


RYY(T,J)E^C(I,J)  -  RYX(I,J)E^C(I,J) 


8 


inC 

where  E  (I,J)  is  the  analytically  specified  E  field  evaluated  at  the  E 

inc  ^  inc  ^ 

mesh  points  and  E  (I,J)  is  the  analytically  specified  E  field  evaluated 
yx  g  g  x 

at  the  Ex  mesh  points.  Finally,  Exy(I,J)  and  Ey^I.J)  are,  in  analogy  with 
eq.  (30), 


1  I 


Exy(I,J)  -  (E^(I,J)  +  E^y(I+l(J)  +  E^(I,J  1)  +  E^y(I+l,J-l))/4  (33) 


E^d.J)  -  (Exxd,J)  +  E^x(I,J+l)  +  Exx(I-l,J)  +  Exx(I-l , J+l) )/4  (3V) 


It  is  also  necessary  to  evaluate  both  components  of  [J  1  as  given  in 

eq.  (9)  at  both  E  and  E  points  in  each  cell.  Using  the  same  convention  as 

x  y 
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x  jriC  'S-il:,'  <zif.  y&S&'i'!*  >,■•&**'  ..*v.*h%V  '*&m^U**r*<*  ^  .<  .  4-«^i>*_, 


17«C,K75ff.  »  V  V  'S'*  KTV^-!  -1 


T  T  T  T 

above,  we  let  J  be  J  evaluated  at  the  E  points  and  J  be  J  evaluated 

xy  y  x  m  yx  x 

at  the  Ey  points.  Then  equation  (21),  where  [J1]  is  actually  required,  uses 
[J^j  in  the  form 


[oofV1]  -  -  [ES]  +  -  Kf  V  x  Hscat] 


S  T 

We  have  already  described  how  to  find  the  [E  ]  contribution  to  [J  ] .  It  is 

easy  to  evaluate  [o0]  ^[Jj]  because  is  a  prescribed  analytic  forcing  term 

which  can  readily  be  evaluated  at  either  the  E  or  the  E  points.  The 

scat  ^  ^ 

troublesome  term  is  [V  x  H  ] .  This  will  have  both  an  x  and  a  y  com¬ 
ponent,  each  of  which  must  be  evaluated  at  the  E  and  the  E  points. 

x  y 

scat 

Let  us  designate  (V  x  H  )  as  the  x-component  of  this  term 
evaluated  at  the  E^  points: 


rv  v  HscaS  _  fld.U.J.J  -..hdU.J- 
K  M  ;xx  Y(J)  -  Y(J-l) 


The  y- component  evaluated  at  an  Ex  point  is 


(V  x  Hscat) 


2L 


Hzq+i.j-i,n  -  (Hzq.jn  ±  Hza.j-n 

2(X(I+1)  -  X(I)) 


(hz(i.j)  +  Hzd.j-m  -  (Hza-i.j')  ±  Hza-i.j-m' 

2(X(I)  -  X(I-l)) 


The  y  component  evaluated  at  an  E^  point  is 


scat.  _  HZq.J)  -  HZ(I-1.J^ 

v  a  'yy  "  X(I)  -  x(i-i) 
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m 


and  the  x  component  evaluated  at  an  E  point  is 

y 


yt_  „scat- 
(V  x  H  ; 

'y x 


i 


a 


irrrizd.j+n  +  Hza-i.j+m  -  (hz(i.j^  -♦-  Hzn-i.jm 
2 [  2(Y(J+1)  -  Y(J) ) 


(HZd.J)  +  HZ(I-l.J))  -  (HZd .  J-l)  -  HZ(I-l.J-l))] 
2(Y(J)  -  Y(J-l))  I 


(39) 


Consequently,  for  example,  eq.  (21)  for  advancing  EX(I,J)  in  all  its 
glory,  becomes 


EX(I,J)n+1/2  -  QXX(I,J)EX(I,J)n'1/2  +  QXY(I,J)EY(I,J)"'1/2 
+  (1  -  QXX(I,j))Ejx(I,J)n  -  QXY(I,J)E^y(I,J)n 
-  SXX(I,J)(Jf(I,J)^x  -  (7  X  HSCat)^) 

*  SXf(I,J)(J-(I,J)"  -  (V  X  aSCat)"  )  (40) 

i.  *-y  Ay 


where  Q XX  and  QXY  are  defined  after  eq.  (30),  EY(I,J)X  is  defined  by  eq. 

(30),  ES  (I,J)  is  defined  by  eq.  (31),  ES  fl.J)  is  defined  by  eq.  (33),  SXX 
xx  xy 

and  SXY  are  defined  after  eq.  (30),  Jf(I,J)  and  Jzr(I,J)  are  the  forcing 

enah  XX  XV 

currents,  (7  x  Hs  c)  is  defined  by  eq.  (36),  and  (7  x  Hs  u)  is  defined 

xx  xy 

by  eq.  (37) . 
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The  scalar  equations  for  advancing  HZ(I,J),  eq.  (20),  is  much  easier  to 
implement  than  the  matrix  equation  for  advancing  [S(I,J)J.  We  now  need  to 
define 

XMUZ(I,J)  as  the  bulk  permeability  at  cell  (I,J),  of  eq.  (4), 

SGMZ(I,J)  as  the  bulk  magnetic  conductivity  at  cell  (I,J),  a0  of  eq. 

(4), 


QMZZ(X.J)  -  s'1™”-')'1  •  SGMZ(I.J)  ■  it 


(41) 


SMZZ(I , J)  -  (1  -  QMZZ(I,J))/SGMZ(I,J)  (42) 

TAUMZZ(I.J)  -  a**1(pa>  -  ^0)  (43) 

evaluated  at  the  center  of  cell  (T,J),  and 

RKZZ(I,J)  -  -  o*)  (44) 

also  evaluated  at  the  center  of  cell  (I,J> 

The  murderously  complicated  interpolations  involved  in  advancing  E  do  not 

occur  in  advancing  H  partly  because  H  is  the  only  component  of  H  present, 

z  z 

and  partly  because  H  is  evaluated  at  the  center,  not  on  an  edge  of  the 
cell. 

g 

From  eq  (10),  H  (I,J)  is  then,  with  frequency  dependence  still  not 

zz 

included, 

HS  (I,J)  -  -  TAUMZZ(I,J)Hinc(I,J)  -  RMZZ(I,J)H^C(I,J)  (45) 

ZZ  ZZ  Z4 

5  tic  inc 

where  H^z  (I,J)  is  the  analytically  specified  Hz  field  evaluated  at  the 

cell  centers. 

One  also  need  only  evaluate  J  (I,J)  of  eq.  (8)  at  the  cell  centers. 
*T  z 

Equation  (20),  where  J  (1,0)  actually  appears,  uses  the  form 

2  15 


<4)‘VTa,J)  -  -  h®z(i,j)  + 


+ 


<4>‘V  x  Iscat) 


zz 


(46) 


Here ,  we  already  have 
current  density  (which 
(V  x  Sscat)zz  is  just 


S  * 

found  H  (I,J),  J-(I,J)  is  a  prescribed  magnetic 
ZZ  jt  zz 

would  be  zero  on  any  physically  real  problem) ,  and 


(V  X  Iscat> 


ZZ 


EY(I+1, J)  -  EY(I,J) 
X0(I+1)  -  X0 (I) 


EX(I , J+l)  -  EX(I,J) 
Y0(J+1)  -  Y0(J) 


(47) 


Thus,  eq.  (20)  for  advancing  HZ(I,J)  becomes 


HZ(I,J)n+1  -  QMZZ(l,J)HZ(I,J)n  +  (1  -  QMZZ(I,J))HS  (I,.I)n+1/2 


-  SMZZ(I,J)(jJ(I,J)^z1/2  +  (V  x  ESCat)ri+l/2)  (48) 
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Introduction  of  Frequency  Dependence 


Let  us  assume  the  frequency -dependent  term  in  eq.  (19)  has  a  kernel  of 
the  form 


M 


K</0 


-I 


a  e 

— m 


-V 


m-l 


(49) 


where  has  ue  units  of  conductivity.  This  assumption  is  equivalent  to 
expanding  the  frequency- dependence  of  the  material's  electrical  propertie.s 
in  a  Prony  series  under  the  constraint  that  all  the  poles  be  on  the  real 
axis.  (Appendix  1  indicates  how  one  may  relax  the  real-poles-only 
constraint. ) 


Equation  (19)  then  becomes 


£ 


3ESCaC 

~~aT~  +  2° 


Escat 


-m 


(50) 


where 


jscat<t) 

-nn 


-P  t  ft  3ESCat(t' )  B  t' 
'ml  _  1  m 

e  J - 


(51) 


T 

In  eq.  (50),  J  is  still  giver,  by  eq..  (35),  but  with  the  understanding  that 

G  S 

E  has  the  frequeney-depe-idont  term  restored.  In  other  words,  E  is  now 

represented  by  eq.  (11),  not  eqs.  (31)  and  (32). 


scat 

The  ^  of  eq.  (50)  are  not  clearly  Identifiable  either  as  conduction 

or  as  displacement  currents.  We  shall  ~oin  the  name  "Prony  currents”  for 

s  c  s  ^ 

them.  Equation  (50)  for  “(t)  is  i^..cb  easier  t'  recognize  if  we  dif¬ 

ferentiate  it  once;  its  homogeneous  solution  is  just  a  decaying  exponential; 
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(52) 


„ Tscat 
dJ 

-=? -  +  B  JSCat 

ot  'm-m 


SE 


scat 


at 


In  the  2D  TM  case,  the  components  of  eqs.  (50)  and  (52)  then  comprise 
2 (M+X )  coupled  first-order  differential  equations  for  Escat  and  Jacat . 
Ideally,  these  equations  should  all  be  advanced  from  (n- 1/2) At  to  (n+l/2)At 
simultaneously  each  cycle.  A  technique  for  doing  this  is  also  described  in 
Appendix  1. 

However,  the  present  code  actually  implements  a  slightly  less  accurate 
sc  at 

algorithm  where  E  alone  is  advancea  first  in  each  cycle,  and  then  the 

scat 

J  are  advanced  separately.  Finally,  a  correction  is  made  to  the  ad- 

^  scat  scat 

vanced  E  '  to  account  for  the  effects  of  the  J 
—  ~m 

At  this  point,  it  is  most  instructive  to  go  back  to  eq.  (7)  and  perform 
a  rearrangement; 


„  dEsM  „ 

7  X  jjscaf-  .  e  .  -^r-  2„  .  £aaa<=  +  .  I«0)  , 


at 


at 


+  (£o  -  1%)  *  Ilnc  +  J  KCt-t') 

-no 


r 


c)ET(t') 

at'  +  h 


(53) 


The  inhomogeneous  part  of  this  equation  can  be  written 


[JT]  +  [JP]  - 


-  [*oHES]  +  [JP]  +  [Jf] 


[V  x  HSCSt]  +  [JP] 


(54'; 


where  [ES]  and  [JP]  are 


[ES]  -  -  [r][EinC]  -  [R][EinC] 


(55) 
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[JP]  -  JC  [K(t-t.')][£T(t')]dt' 


with  [r]  and  [R]  still  respectively  given  by  eqs.  (15)  and  (16). 


Substituting  eqs.  (54)  -  (56)  back  into  (53)  then  results  in 


-  -scat  .T  .P 

A.  •  ~sT~  +  2o  '  E  -1 


s  ^ 

K  P$ 

S  £ 


P  f? 

5  4? 


This  equation  is  just  (19)  with  the  frequency -dependent  term  transferred  to 

P 

the  right  and  represented  as  J  .  It  is  advanced  according  to  eq  (21): 


n+1/2  “t€oJ  [°o]dt  n-1/2 

[E(I,J)]n+i/Z  -  e  [E(I,J)]n  L/Z 


\  1 


*teoo^  t»0]At  -IT  P  n 

-  ([I]  -  e  )[a0]  1[Ji(I,J)  +  JP(I.J)]n  (58) 


The  problem  with  direct  application  of  this  procedure  is  that  we  do  not 

P  •  S  Celt 

know  the  portion  of  J  associated  with  E  at  nAt  until  we  have  advanced 
the  Prony  currents,  and  we  cannot,  strictly  speaking,  do  that  until  we  have 
advanced  Escat. 


As  mentioned  previously,  the  code  does  not  presently  utilize  the  proce- 

scat 

dure  described  in  Appendix  1  for  simultaneous  advancement  of  E  and 

scat  scat 

J  Rather,  we  first  find  an  intermediate  value  of  E  obtained  with 

“m  _ 

effects  of  the  Prony  currents  omitted: 


int  "t€co^  MAt  n-1/2 

[E(I , J) ] int  -  e  [E(I,J)]a  L/Z 
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tt73t.rJLTtT;P&*X?S.'  N  *  ^araHTOV^  -i^TCTTWTOroTataffTCira  xxjm**^* 


I 


8 


« 


j 

! 


I 

N 

K 

S 

I 


1 

i 


s 

ft 


-[*„,]  [ffoJAt  -IT  n 

-  (EX]  -  e  )[a0)  1[JA(I,J)]xl  (59) 


The  procedure  for  obtaining  this  intermediate  value  is  identical  to  the 
procedure  described  in  the  previous  section  for  advancement  through  a  total 
cycle  in  the  absence  of  frequency -dependent  effects.  Its  implementation  in 
the  code  is  also  identical  to  what  was  described  in  the  previous  section. 

P 

Let  us  next  turn  to  the  advancement  of  the  total  Prony  current  J  as 
given  by  eqs.  (49)  and  (56): 


M 


m-1 


(60) 


where 


JT(t) 

—ra 


-0mt  ft  3ET(t' ) 

j  .___ 


A-t' 

e  dt' 


(61) 


Equation  (61),  like  (51),  is  made  more  recognizable  by  differentiating  with 
respect  to  time: 


T 

y9  J 

'm-m 


9ET  3(Elnc  +  EScat) 
3t  ”  3t 


(62) 


The  equation  for  exponential -difference  advancement  of  this  result  is 


JT(I,J)n+1/2  -  e 
-m 


-fl  At  rr  ,  n 

“  J1(T,J)Tl'1/2 

"nr  ’ 
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'  W  If  v  •"*.**?.  fBICTBtWMfWW'rTTT'FT’rf  VI 


(63) 


-A.A* 


+  (1  -  e  m  )<iinC(I,J)  +  iSCat(I,J))n/i3„ 


p 

I 

b 


I® 


The  incident  field  £^nc  is  a  specified  analytic  function.  At  present,  the 
code  evaluates  fescat(I,J)n  as 


|5S 


W. 


iscat(i,j)n 


S(I,J)int  -  I(l,j)n‘1/2 


At 


(64) 


8 

m 


§ 

© 


SR 

'*1 


I 

W 


I 


S 


H 


Equations  (63)  and  (64)  nay  be  combined  to  give  an  expression  for 
J^(I,J)n+^/2  in  terms  of  known  quantities.  We  can  ther  determine  JP(I , J)n 
as 


ra.j) 


n 


M 


}  [am(I,J)][J^(I,J)n+1/2 

m=l 


+  J*(I,J)n'1/2]/2 


(65) 


Subtraction  of  eq.  (59)  from  eq.  (58)  then  permits  us  to  advance 
[E(I , J) ]  from  its  Intermediate  value  to  itc  value  at  the  new  time  step, 
(n+l/2)At: 

[E(I,J)]n+1/2  -  tE(I,J)]int 


-  ([I] 


-uaor1[<703At 

e  )[o0] 


1[JP(I,J)]n 


(60) 


Frequency  dependence  in  the  magnetic  properties  of  the  material  can  be 
treated  in  an  exactly  dual  manner  to  what  we  have  just  described  for  the 
electrical  properties. 
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At  present,  we  have  not  coded  up  magnetic  frequency  dependence,  nor 
have  we  combined  frequency  dependence  with  off-diagonal  type  anisotropy. 
Both  these  generalizations  would  be  perfectly  obvious  extensions  of  what  has 
been  done,  but  we  cannot  conceive  a  canon ’cal  problem  we  could  check  the 
results  against. 

Not  mixing  off-diagonal  anisotropy  and  frequency  dependence  means  we 
only  treat  diagonal  affl  tensors.  If  eq.  (65)  is  substituted  in  eq.  (66),  we 
obtain 


[E(I,J)]n+1/2  -  [E(I , J) ] int  -  ([I]  -  e  “ 

•  ^  [SAm(I,J)][J^(I,J)n+1/2  +  J^(I,J)n_1/2 


m=l 


[o0]At 

) 


] 


(67) 


I 


§5 

& 


where 


[SAm(I,J)]  -  [a0(I,J)]'1[Am(I,J)] 


(68) 


KS 


$ 


9 


•o 

S 


In  the  actual  code,  two  arrays  are  used  to  describe  the  material  ef¬ 
fects  of  each  term  in  the  Prony  series.  Let  represent  the  xx 

element  of  the  inverse  of  the  matrix  described  by  eqs.  (22)  -  (25)  at  cell 
(I,J).  Let  cr0(I,J)yy  similarly  represent  the  yy  element  of  the  inverse  of 
the  matrix  described  by  eqs,  (26)  -  (29)  at  cell  (I,J). 

Moreover,  let  AAMX(I,J)  represent  the  xx  element  of  the  mth  Prony 

tensor  of  the  bulk  material  at  cell  (I,J),  and  let  AAMY(I,J)  represent  the 

corresponding  yy  element.  Then  the  xx  element  of  [SA  (I,J)]  which  actually 

X  m 

relates  the  x  component  of  J  (I,J)  to  the  x  component  of  E(I,J)  is  called 
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SAMX(I,J)  -  ff0(I.J)^x  (AAMX(I.J-l)  +  AAMX(I,J))/2 


Similarly,  the  yy  element  of  [SA  (I,J)]  which  actually  relates  the  y  com- 

X  ^ 

ponent  to  the  y  component  of  E(I,J)  is  called 


SAMY(I.J)  -  cr0(I,J)^  (AAMY(I-l.J)  +  AAMY(I , J) )/2 


In  keeping  with  our  simplification  of  not  mixing  off-diagonal  anisotropy 
with  frequency  dependence,  we  ignore  any  possible  off-diagonal  nonzero 
values  in  [SA  (I,J)J. 

It  turns  out  that  only  SAMX(I,J)  and  SAMY(I,J)  need  actually  be  stored. 
That  is,  it  is  not  necessary  to  assign  arrays  for  keeping  ct0(I,J)^, 
a0(I,J)^,  AAMX(I,J)  and  AAMY(I , J)  . 

Consequently,  the  actual  equation  used  in  the  code  for  implementing  the 
x- component  of  the  Prony  correction  is 


EX(I,J)n+1/2  -  EX(I,J)int 


-  (1  -  QXX(I , J) )  )  SAMX(I,J)  XJMSX(I , J) 


where 


XJMX(I,J)n  -  [J^(I,J)"+1/2  +  j3(I,J)?'1/2]/2 
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Similarly,  the  actual  equation  used  for  implementing  the  y-component  of  the 
Prony  correction  is 


§§ 


f-  ^ 

1  & 

l 


EY(I,J)n+1/2  -  EY(I,J)int 


M 

-  (1  -  QYV(I.J))  ^  SAMY(I.J)  XJMY(I,J)n 


m-1 


where 


XJMY(I,J)n  -  [JT(I,J)n+1/2  +  jT(I,J>n-1''2]/2 
m  y  m  ‘  y 


I 


!§ 


m 

Si 


i 


sa 


m 
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Transformation  from  the  Near-field  Time  Domain  to  the 
Far -field  Frequency  Domain 


The  foregoing  work  described  the  determination  of  the  scatterer's 
electromagnetic  response  and  associated  near  fields.  We  are  actually  inter¬ 
ested  in  the  RCS,  which  is  a  far-field  quantity.  Now  we  shall  describe  how 
the  code  extracts  **he  RCS  from  the  near- field  results.  In  this  process,  we 
also  transform  from  time  domain  to  frequency  domain. 

Any  electromagnetic  field  can  be  expressed  in  terms  of  an  electric  and 

•fa 

a  magnetic  vector  potential,  A  and  A  .  These  vector  potentials  (in  the 
frequency  domain)  obey  the  inhomogeneous  wave  equations. 


V2A  +  k2A  -  -/i0J  (75) 

V2A*  +  k2&*  -  e0l  (76) 


Here  J.  is  the  fictitious  magnetic  current  density  often  found  useful  in 
manipulating  Maxwell's  equations.  Equations  (75)  and  (76)  can  be  general¬ 
ized  to  apply  tc  any  linear  medium,  although  we  shall  find  their  free-space 
form  adequate  for  our  uses. 

In  2D,  we  define  the  far  field  to  be  the  region  where  all  fields  drop 
-1/2 

off  as  r  '  ;  I.e.,  where  all  the  faster  falling  terms  have  vanished.  We 

can  then  separate  the  electric  and  magnetic  fields  into  two  parts, 


E  - 


E 

~e 


h  F. 


~m 


(77) 


H 


H  +  H 
_e  ~ra 


(78) 
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The  e- subscripted  parts  are  associated  with  the  electric  vector  potential  4, 
and  the  m- subscripted  parts  are  associated  with  the  magnetic  vector  poten- 
tial  A  .  In  particular,  if  we  call  Y0  and  Z0  the  admittance  and  impedance 
of  free  space,  we  can  show  in  the  far  field  that 


He  -  V  X  A//i0  -  icoY0ir  X  4  (79) 

^  -  V  x  &*Ao  -  iwZoj^  x  A*  (80) 

— e  "  iw^t  “  iw(VV  +  izAz>  (81> 

\  -  -  i«4*  -  -  1»U^aJ  +  izA*)  (82) 


Here,  a  t  subscript  indicates  that  only  the  transverse  components  (tf>  and  z) 
are  retained.  Equations  (79)  -  (82)  are  analogous  to  3D  formulas,  and 
depend  on  the  fact  that 


€l 


SB 


i(kr-wt) 

c 

Jr 


(83) 


9 

§5 


is  a  valid  far- field  frequency-domain  2D  solution  of  the  wave  equation  even 
if  the  more  general 


f(t  -  r/c 

Jr 


(84) 


I 


1 

5? 


is  not  a  valid  time-domain  solution.  These  equations  tell  us  that  if  we  can 
•k 

evaluate  A  and  A  ,  we  can  find  the  2D  RCS  without  undue  complication. 
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4. 


In  two  dimensions,  the  Green's  function  for  the  scalar  wave  equation  in 
the  frequency  domain  obeys 


V2G(£|r')  +  k2G(rjr')  -  6(tc- r ' ) 


(85) 


where  r  is  the  scatterer  location  and  r'  is  the  observer  location.  This 
equation  has  solution,  with  R  -  |r'  -  r| , 


G(r|r')  -  -  Jh^OcR) 


(86) 


Thus,  at  least  in  cartesian  coordinates,  A  and  A  become 


friH(1)(kR) 

A(£'  ,w)  -  jj — -  I(£,w)d£ 


(87) 


*  ffiH^1)(kR)  * 

A  (r',w)  -  -  JJ - ^ -  J  (r,w)dr 


(88) 


For  the  far  field  region,  G(r|r')  asymptotically  approaches 


G(r|r') 


e-i3ir/4 

4 


(89) 


This  expression  may  be  further  manipulated  by  letting  r'  replace  R  in  the 
denominator  of  the  radical.  The  phase  term  requires  a  bit  more  care: 


kR  -»  k(r'  -  i^  »  r)  -  kr'  -  (kx  cos <f>'  +  ky  sin <f>') 


(90) 
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where 


i,.,  “  i^cos^'  +  iy 


is  a  unit  vector  pointing  from  the  target  to  a  far-field  observer. 

Using  these  expansions,  we  can  rewrite  the  formula  for  A  in  the  far 
field  as 


A(r' ,«) 


-3ijr/4 

M(>e  /  2  ikr 

4  V  wkr'  e 


'  ik(x  cos<t>'  +  *  sin*r)dr 


A  corresponding  expression  exists  for  4  (£',£*>)•  The  far-field  expression 


for  E  then  becomes 
~e 


Eft(r' ,w)  -  iw  At(r' ,w) 


-3iw/4 


hr  «ikr'  If  i.<I,»)e-ik<x  cos*’  *  y  sin#,)dr 


Similarly,  becomes 


HjjjCl'.w)  -  i«At(r',w) 
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(94) 


i  w£0e 


-3i?r/4 


•/ 


irkr  ■ 


ikr' 


JJ  i*(i,»)«'lk(x  cos*' 


+  y  sin^')d£ 


Analogous  formulas  exist-  f-r  E  and  H  . 

~in  e 

Equations  (93)  and  (94)  are  not  directly  applicable  to  the  output  of 
our  2D  Maxwell  solver  as  these  equations  demand  the  frequency- domain  J  and 
J  ,  while  the.  Maxwell  solver  outputs  the  time-domain  currents. 


Let  us  say  we 


«t->t  E  (r'w  )  where  r'  points  in  one  of  N  discrete 
-e'—p’  q  ~p  p 

angles  of  interest  and  is  one  of  discrete  frequencies  of  interest.  We 
can  then  write 


E  (r\u  )  - 
~e~p  q' 


8jt 


Virkr' 

q  p 


(95) 


where  k'  is  the  wavenumber  pointing  towards  the  observer  at  location  r'  and 
~pq  HP 

frequency  and  where  we  have  interchanged  the  order  of  time  and  space 

integration  after  replacing  J.t(r,w)  with  its  inverse  Fourier  representation. 
Analogously,  H^r^.w^)  becomes 


H  (r'  ,o>  ) 
“m  — p  q 


/XT  .Vi 


8?r 


'  rck'r' 

q  P 


J7  v  dt  JJ  4<i, 


-ik'  •  r 
t)e  ~ ^  dr 


(96) 
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Equations  (95)  and  (96),  and  the  two  companion  equations  for  E^  and  Hg 
are  in  a  form  which  is  compatible  with  our  time -domain  Maxwell  solver. 
Taking  the  time  integration  outside  the  space  integration  is  vitally  impor¬ 
tant  to  the  efficiency  of  our  algorithm.  Were  this  not  done,  J  (r,t)  and 
J^(r ,  t)  would  hace  to  be  Fourier  transformed  at  every  point  before  being 
integrated  over  space.  The  form  of  eqs.  (95)  and  (96)  replaces  this  enor¬ 
mous  computation  with  a  single  Fourier  transform  on  the  result  of  the  space 
i ntegral . 

SCclt 

If  we  let  J  represent  the  total  current  (conduction,  displacement 
and  Prony)  associated  with  the  scattered  electromagnetic  field  (see  eq. 
(53)). 


f 


..scat 

J  —  e  • 


3ESCat  5EinC 

+  £o  •  ISCat  +  (U  '  1* o) 


at 


at 


3ET(t' ) 


.  OCj  { C.  ) 

+  (£o  -  1%)  •  IlnC  +  J  K(t-t')  .  — — —  dt' 


(97) 


SC3.tl  • 

and  If  we  substitute  J  for  J  in  cq.  (95),  E  of  eq.  (95)  becomes  the 

.  C  C  6 

scat 

scattered  field  E  of  the  first  section  unless  magnetic  materials  are 

s  cat 

present.  That  is,  J  integrated  over  the  scatterer  cross-section  accord- 

SC3C 

ing  to  eq.  (95)  gives  E  in  the  absence  of  magnetic  materials. 


pa 


RCS  (*',«  )  =  2?r 

p  q 


..scat.  ,  . 

E  (r' ,w  ) 
de _ ±P-l_a_ 


E1(r' ,w  ) 
-  v-p  q' 


(98) 


In  this  convention,  A  and  E  are  zero. 

—  ~m 

However,  it  is  possible  to  replace  the  area  integral  of  eqs.  (95)  and 
(96)  with  a  contour  integral  by  means  of  Huygens  principle.  In  particular, 
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let  S  be  a  closed  contour  which  completely  surrounds  the  target.  For  in¬ 
stance,  let  S  be  a  rectangle  defined  by  x  -  x*^  and  y  -  ylty2. 

(At  this  point,  it  may  be  well  to  digress  a  moment  and  present  a  map  of 
the  problem  space  we  are  using.  Let  us  refer  to  Figure  2.  The  entire 
problem  space  is  described  by  1  i  I  S  NX,  1  <  J  ^  NY  or  X0(l)  <  x  £  X0(NX), 
Y0(l)  <,  y  :£  Y0(NY),  where  NX  and  NY  are  typically  around  250.  The  working 
volume  actually  occupied  by  the  scatterer  is  (NXB  -  1)  by  (NYB  1)  cells 
centered  in  the  overall  problem  space  where  NXB  and  NYB  are  typically  on 
the  order  of  75.  Thus,  the  working  volume  is  separated  from  the  problem 
space  boundary  by  (1/2) 'NX  •  NXB)  cells  along  x  and  (1/2) (NY  -  NYB)  cells 
along  y.  This  separation,  which  is  on  the  order  of  80  culls,  is  necessary 
to  decouple  the  outer  boundary  from  the  reactive  fields  of  the  scatterer. 
The  Huygens  surface  S  is  normally  placed  me  cell  outside  the  actual  working 
volume,  at  I  -  ILOW  or  IHIGH  and  J  -  JLOW  or  JHIGH .  Then  xx  is  Xo(IL0W), 
etc. ) 


sc&ti  S(at 

Let  5  and  H  be  the  scattered  fields  which  our  Maxwell  solver 
predicts  will  exist  on  S  due  to  the  time -domain  illumination.  Let  n  be  an 
outward-pointing  unit  normal  on  S.  If  we  then  remove  the  scatterers  and  its 
currents,  but  let  an  electric  surface  current 


K  -  n  X  H 


scat 


(99) 


and  a  magnetic  surface  current 


* 


K 


n  X  E 


scat 


(100) 


flow  oh  S, 
This  means 


the  scattered  electromagnetic  field  will  be  replicated  outside  S. 
the  area  integral  of  eq,  (95)  may  be  replaced  by 
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Map  of  overall  problem  space  showing  location 
working  volume  and  Huygens'  surface. 


>'5£l  d^^Jrtlk 


-ikV  •  r. 


e  (a  X  fisoat(ii,t))te  -M  -1  a.t  =  j£«<t> 


<101) 


where  the  summation  over  i  represents  integration  over  the  finite  difference 

eraf*  SCa.t 

cell  edges  which  lie  on  S  and  H  (r^,t)  is  II  evaluated  at  the  middle 

of  As^.  This  summation  is  represented  by  I  because  it  has  the  dimension  of 

amperes.  In  the  actual  code,  it  goes  by  the  name  XIEST^  if  one  is  computing 

a  monostatic  RCS  or  XIESTB  if  one  is  computing  a  bistatic  RCS.  Figure  1 
scat  * 

indicates  that  H  is  actually  evaluated  at  the  cell  centers,  not  on  the 

scat 

cell  edges.  Thus,  an  interpolation  is  necessity  to  obtain  H  (r^,t).  For 
instance,  on  the  y  -  Y0(JLOW)  portion  of  S,  we  have 


(n  x  Hscac(ri,t))t  -  -  i^  sin  (HZ(I,JL0W)  +  HZ(I,JL0W  -  l))/2  (102) 


Analogously,  the  area  integral  present  in  eq.  (96)  for  H^(r^.ui^)  can  be 
written 


S  (-n  X  iscat(Ei,t))te  1 

i 


As,  -  Iscat(t)* 
i  t 


(103) 


(Remember  the  t  subscript  on  I,  I*,  or  anything  else  implies  evaluation  with 

the  r  component  omitted.)  In  the  code,  this  variable  is  called  XIMST  or 

scat  * 

XIMSTB  .  Since  E  is  evaluated  at  the  middle  of  the  cell  edges,  no 

F  scat 

interpolation  is  necessary  to  obtain  E  (r^,t).  On  the  y  -  Y0(JLOW) 

portion  of  S,  the  equation  analogous  to  (102)  is 


(-n  X  ES  (r. ,t)).  -  -  i  EX(I.JLOW) 
1  t  z 


(104) 


Equations  (95)  and  (101)  indicate  the  z-component  of  gscat  £s  as. 

sociated  with  the  ^'-component  of  E  (£,«  ).  The  ^-component  of  Jg  in  eq. 

e  q 

(103)  has  the  same  association,  as  the  cross  with  i  in  eq.  (80)  proves. 
These  scattered  fields  are  the  TM  solution. 


SCdtl 

Analogously,  the  ^-component  of  K  in  eq.  (101)  and  the  z-component 

SCclt 

of  E  in  eq.  (103)  relate  to  the  (decoupled)  TE  solution,  which  we  are 

not  treating  in  detail  at  this  time. 


Let  us  use 


to  denote  this  "current"  evaluated  from  the 


iscat  n 

“Pq  t  *n+l/2 

finite -difference  code  at  nAt.  Let  us  analogously  denote  1^  (t)t  '  . 

As  the  finite -difference  calculation  progresses,  we  can  then  keep  running 

summations  of  each  scattering  direction  and  each  frequency  w^, 


0scat 

^pq 


n 

2 

m-1 


jScat 

n?q 


ju  mAt 
e  q  At 


(105) 


Qscat 


(t)*n+l/2 


n 

2 

m-1 


Tscat,  »*m+l/2 

npq  ^  't 


jw  (m+l/2)At 
e  q  At 


(106) 


When  the  time-domain  finite  difference  calculation  is  complete,  these  Q's 

will  then  respectively  contain  quantities  which  are  directly  proportional  to 

the  electric  and  magnetic  contribution  to  the  RCS  in  direction  r'  or  at  w  . 

p  q 

Note  that  it  is  only  necessary  to  back  store  2(N^  +  N^)  complex  quantities 
during  the  time-domain  finite-difference  calculation  in  order  to  preserve 
all  the  information  necessary  to  generate  a  monostatic  RCS  as  a  function  of 
o)  and  a  bistatic  RCS  as  a  function  of  scattering  direction. 


The  symbol  Q  is  used  in  eq.  (105)  because  it  represents  a  quantity  with 
units  of  coulombs.  In  the  code,  it  is  written  SIEST^  if  one  is  computing  a 
monostatic  RCS  or  SIESTBp  if  one  is  computing  a  bistatic  RCS.  The  scattered 
electric  field  associated  with  Q^qat(t)t  is 
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,  -3iw/4 


K  _ uae  '  r~r  ■  r  ik  r' 

Escatfr'  «  )  -  -  — ^ -  / -2 -  e  q  P  0scatf»)® 

£e  (Ep  q;  8*  v  *k  e  Sq  wt 


(107) 


The  code  tracks  the  ^'-component  of  this  quantity  as  EPES^  or  EPESB^,  when 

Jr'  is  normalized  out.  Analogously,  the  scattered  magnetic  field  associated 
P  _scat/4..*  . 

Wlth  2pq  <t)t  iS 


Hscat(r,  j 
-m  v— p  q 


,  -3iw/4 

iu>a£oB  pi 
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\  *k  r' 

q  P 


ik  r' 

.  q  P  n3cat(«)*e0 
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(108) 


In  the  code,  the  z- component  of  this  quantity  is  HZMS^  or  HZMSBp  when  is 
normalized  out.  The  scattered  electric  field  associated  with  the  magnetic 
current  is  obtained  by  combining  eqs.  (80)  and  (108), 


—SC^t*  rt  \ 

E  (r',u  )  - 
~m  ^p  q 


.  -3i*/4 

^a6  Hi 


8jtc 


*k  r' 
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,  q  P 


I^p  x  Q®cat(»)*“  (109) 


In  the  code,  the  ^'-component  of  this  quantity,  less  the  J r^  factor,  is 
EPMS  or  EPMSB  . 

q  p 

The  radar  cross  section  then  becomes 


RCS(^' ,w  ) 

p  q 


2n 


(Escac(r',w  )  +  Escac(r' ,w  ))/r' 
v~e  v~p’  q '  ~p  a 


E1(r ' ,u  ) 

-  v-p’  q/ 


(110) 


In  the  code,  SCATXE  or  SCATXEB  is  the  E  contribution  to  the  ratio 

q  P  e  scat 

inside  the  absolute  value  signs,  and  SCATXM  or  SCA1XMB  is  the  E 

°  q  pm 

contributions . 
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-  Y0i2  J0(k0qr)  +  ^  2Jn(k0qr)cos  n<£  +  ^  2iJn(k0qr)sin  n<£  (3) 

n-2 , 2  n-1 , 2 


In  the  future,  it  will  be  useful  to  designate  the  coefficients  of  these 


harmonics  as  a^nC; 


n  >  0,  even 


n  odd 


i-Tic  inc  inc 

Since  V  X  H  -  -  juqe0E  ,  the  cylindrical  harmonic  expansion  for  E 

becomes 


r-  CO  CO 

iY  i 

_inc .  .  0  ~r  \  inc. . inc. 

E  (r.u  )  —  -  -  —  )  a  J  (kn  r)n  sin  nd>  -  )  a.  J 

-  q'  w  e0  r  /  n  nv  Oq  '  T  [_  n  n 

q  [  [n-0,2  n-1,2 


(kQqr)n  cos  n$ 


\  a^ncJ'(kn  r)cos  n<j>  +  \  a^ncJ' 
(_  n  n  Oq  Y  [_  n  n 

n-0 , 2  n-1 , 2 


(h0qr)sin  n<f> 


The  innermost  material  will  include  the  cylinder  axis.  Thus,  only 
Bessel  functions  of  the  first  kind  are  permitted  in  the  solution  there: 


HKr.cOq) 


CO  CO  -» 

Yn  i  )  a1  J  (k.  r)cos  n <f>  +  )  a1  J  (k,  r)sin  n<f> 

lq~z  j  nq  nv  lq  ‘  Y  /  nq  n'  lq  r 

n=0 , 2  n-1 , 2 
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)  a1  J'(k,  r)cos  n<£  + 
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n*0 , 2 
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where  Y,  is  the  admittance  of  medium  1 

iq 


l*i  +  JgiA\ 

'  Mi  +  jo'i/Wq 


(8) 


and  k,  is  the  wavenumber  of  medium  1  at  w  , 
lq  q’ 


(9) 


The  N-l  concentric  shells  will  permit  solutions  of  both  kinds.  Thus, 
in  region  i,  1  <  i  <  N,  of  eqs.  (6)  -  (9)  becomes  replaced  by 


a1  J  (k,  r)  -»  a"  J  (k,  r)  +  Y  (k.  r) 
nq  r  lq  nq  n  iq  nq  n'  iq 


(10) 


Finally,  in  free  space  outside  the  cylinder,  the  first  Hankel  function 
is  the  only  permitted  solution  for  the  scattered  field.  Thus,  in  this 

inc 

region,  J^CkQ^r)  is  replaced  by 


me  ■»  \ 

a  J  (kA  r) 
n  n  Oq 


aScatH<l)(I(  r) 
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(11) 
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in  eqs.  (3)  -  (5) . 

The  boundary  conditions  at  each  interface  are  that  eEr,  and  Hz  be 
continuous.  It  turns  out  that  the  first  and  third  of  these  conditions  are 
equivalent.  Thus,  matching  of  coefficients  at  the  innermost  interface  lead: 
to 


J?' 

# 


a1  Yt  J  (k.  a.)  -  a2  Y„  J  (k„  a.)  -  b2  Y„  Y  <k„  a! )  -  0 
nq  lq  n^  Iq  V  nq  2q  nv  2c  V  nq  2q  nv  2q  1' 
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a1  Jf (k-  a,)  -  a2  J'(k„  a,)  -  b2  Y' (k0  a,) 
nq  nv  lq  V  nq  nv  2q  V  nq  n'  2q  V 


(12) 


Matching  of  coefficients  at  any  other  interface  except  the  outer  bound¬ 
ary  of  the  cylinder  gives 
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a1  J ' (k.  a.  .)  -  b1  Y'(k.  a.  .)  -  0 
nq  n  iq  l-l  nq  n  iq  i-l 


(13) 


Finally,  the  boundary  condition  at  the  outermost  surface  is 
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+  -  CtMn1)'(kOqa»)  ‘  4“Ji(k0,aS> 


(14) 


For  each  azimuthal  harmonic  and  each  frequency,  eqs.  (12)  -  (14)  com¬ 
prise  a  set  of  2N  linear  equations  in  2N  unknowns.  The  associated  matrix  is 
five-banded,  and  extremely  easy  to  solve  by  Gaussian  elimination.  (The  main 
diagonal  and  first  diagonal  off  each  side  of  the  main  is  full.  The  second 
diagonal  off  each  side  of  the  main  is  half  zeros.) 


The  quantities  of  interest  in  RCS  evaluation  are 
dimensional  bistatic  RCS  is  defined  by 
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The  scattered  electric  field  (neglecting  the  reactive  radial  component)  is 
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Here,  use  is  made  of  the  identity  ^QqY0  *=  w^e0. 


For  large  arguments,  Hankel  functions  have  the  asymptotic  limit 
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- *> 


(17) 
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Substitution  of 
scat 


the  a~ 


nq 


eqs.  (16)  and  (17)  in  (15)  then  yields  the  RCS  in  terms  of 


RCS($* ,«  ) 
P  q' 


4_ 

k. 


n-0,2 


scat.  ..n  ,, 

a  ( -  a )  cos  no 
nq  *p 


n-1,2 


scat  #  . >n  >  . , 

anq  (‘1)  sin  % 


(18) 


A  total  of  six  code-code  comparisons  have  been  made  between  the 

cylindrical -harmonic  frequency- domain  algorithm  and  out  TDFD  code.  First  we 

examined  a  solid  dielectric  circular  cylinder  of  e  -  2.  The  cylinder  was 

inc 

given  a  radius  a  -  .5  m.  We  used  an  angle  of  incidence  of  <j>  -  +  45°  from 

the  x-axis.  Meshing  was  square,  with  each  cell  .04  m  on  a  side.  Figure  3 
illustrates  the  TDFD  cross-sectional  model  for  the  cylinder. 

Figure  4  is  a  linear  comparison  of  the  monostatic  result,  and  Figure  5 
is  a  dB  comparison.  Figure  6  is  a  linear  comparison  of  the  bistatic  result 
of  250  MHz,  and  Figure  7  is  a  dB  comparison.  Figures  8  and  9  are  bistatic 
comparisons  at  500  MHz. 

The  second  comparison  was  an  identical  cylinder  except  that  e r  was 
increased  to  -  9.  Figures  10  -  15  compare  the  same  data  as  Figures  4-9 
did  for  the  first  example. 

Figure  16  is  a  four-way  dB  comparison  of  the  monostatic  RCS  computed  by 
both  techniques  for  both  values  of  e  . 

Good  agreement  occurs  for  all  the  er  -  2  comparisons,  but  discrepancies 
occur  above  300  MHz  when  e  —  9.  For  instance,  agreement  is  not  good  in  the 
bistatic  result  at  -  9  when  frequency  is  500  MHz  (Fig.  15).  If  e  *»  9, 
300  MHz  corresponds  to  a  wavelength  in  the  cylinder  of  .33  m,  or  about  16 
cells  per  wavelength.  The  conclusion  that  a  wavelength  must  be  resolved 
into  16  segments  to  get  good  finite-difference  results  is  neither  new  nor 
surprising. 
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DIELECTRIC  ROD  CROSS  SECTION 


Figure  3.  TDFD  model  of  a  dielectric  cylinder  with  . 
radius  and  .02  m  square  cells. 


Monostatic  Radar  Cross  Section 


Monostatic 


Bistatic  RCS  for  a  circular  dielectric  cylinder 
with  «  -  2  at  250  MHz.  Solid  curve  is  the  TDFD 

result,  and  dashed  line  is  the  cylindrical - 
harmonic  result.  Plot  is  based  on  illumination 


Dielectric  Rod  Bistatic  Radar  Cross  Section 


Bistatic  RCS  for  a  circular  dielectric  cylinder 
with  er  -  2  at  500  MHz.  Solid  curve  is  the  TDFD 

result,  and  dashed  line  is  the  cylindrical - 
harmonic  result.  Plot  is  based  on  illumination 


Bistatic  RCS  for  a  circular  dielectric  cylinder 
with  e  -  9  at  250  MHz.  Solid  curve  is  the  TDFD 

result,  and  dashed  line  is  the  cylindrical- 
harmonic  result.  Plot  is  based  on  illumination 


Dielectric  Rod  Bistatic  Radar  Cross  Section 


Bistatic  RCS  for  a  circular  dielectric  cylinder 
with  €  -  9  at  500  MHz.  Solid  curve  is  the  TDFD 

result,  and  dashed  line  is  the  cylindrical- 
harmonic  result.  Plot  is  based  on  illumination 


Dielectric  Rod  Bistatic  Radar  Cross  Section 


Bistatic  RCS  for  a  circular  dielectric  cylinder 
with  e  -  9  at  500  MHz.  Solid  curve  is  the  TDFD 

result,  and  dashed  line  is  the  cylindrical- 
harmonic  result.  Plot  is  based  on  illumination 


RCS  of  a  circular  dielectric  cylinder 


I 
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The  third  comparison  was  for  a  perfectly  magnetically  conducting 
cylinder  of  a  -  .5  m  radius.  Meshing  was  again  square,  but  with  each  cell 
now  .04  m  on  a  side.  Figure  17  illustrates  this  TDFD  model. 

(TM  scattering  off  a  perfectly  magnetically  conducting  cylinder  will 
give  identical  RCS  results  as  the  more  familiar  problem  of  TE  scattering  off 
a  perfectly  electrically  conducting  cylinder.) 

Figure  18  shows  a  linear  comparison  of  the  two  monostatic  RCS  calcula¬ 
tions,  and  Figure  19  gives  a  dB  comparison.  Figure  20  is  a  linear 
comparison  of  the  bistatic  result  at  250  MHz,  and  Figure  21  is  a  dB 
comparison.  Figures  22  and  23  are  bistatic  comparisons  at  500  MHz. 


a 


in 

*  » 1 


i 


I 
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Although  all  these  magnetic- cylinder  results  show  good  agreement  be¬ 
tween  the  two  techniques  an  interesting  minor  difference  does  appear, 
especially  in  the  second  and  fourth  quadrants,  of  Figures  21  and  23.  In 
particular,  the  cylindrical -harmonic  bistatic  RCS  results  are  exactly  sym¬ 
metrical  about  the  main  diagonal  ( <f>  -  45°),  while  the  TDFD  result  is  not. 
We  believe  a  low-grade  bug  is  present  in  the  TDFD  coding,  but  that  this  bug 
was  undetectable  until  bistatic  capability  was  added  to  the  TDFD  code  in  the 
final  days  of  the  effort. 

The  fourth  comparison  was  for  a  perfectly  magnetically  conducting  rod 

of  the  same  radius,  but  coated  by  a  damper  .5  m  thick.  Thus,  the  overall 

radius  of  the  composite  rod  was  now  1  m.  The  damping  material  had 

.3 

properties  «r  *■  1,  /*r  -  1,  o-e  o*/p  -  4  x  10  .  These  values  were 
selected  to  give  a  skin  depth  on  the  order  of  the  damper  thickness  at  fre¬ 
quencies  (50  -  500  MHz)  for  which  calculations  were  run.  Square  cells  .04  m 
on  a  side  were  again  used. 

Figure  24  illustrates  the  TDFD  model  of  the  damped  cylinder..  Figures 
25  -  30  compare  the  same  data  as  Figures  18  -  23  did  for  the  bare  cylinder. 

It  is  important  to  note  that  damping  the  rod  reduces  the  RCS  in  all  but 
the  forward  direction*  where  it  enhances  the  RCS.  This  occur  because  the 
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Figure  17.  TDFD  model  of  a  magnetically  conducting 
cylinder  with  .5  m  radius  and  .04 
cells . 


^J*.^  -  if’i  *ZXi7t^Z  --‘w-  .T.~  ^  .  • 


Figure  18.  Monostatic  RCS  for  a  bare,  magnetically  conduct¬ 
ing  circular  cylinder.  Solid  curve  is  the  TDFD 
result,  and  dashed  line  is  the  cylindrical- 
harmonic  result.  Plot  is  based  on  illumination 


Figure  24.  TDFD  model  of  a  magnetically  conducting  cylinder 
of  .5  m  radius  coated  by  a  damper  .5  m  thick. 
Cells  are  .04  m  square. 


Coated  Magnetic  Rod  Monostatic  Radar  Cross  Section 
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Coated  Magnetic  Rod  Bistatic 


Coated  Magnetic  Rod  Bistatic  Radar  Cross  Section 
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Magnetic  Rod  Bistatic  Radar  Cross  Section 


0  MHz.  Soli' 
dashed  line 


Coated  Magnetic  Rod  Bistatic  Radar  Cross  Section 


30  Bistatic  RCS  for  a  coated,  magnetically  conduct¬ 
ing  circular  cylinder  at  500  MHz.  Solid  curve 
is  the  TDFD  result,  and  dashed  line  is  the 
cylindrical  harmonic  result.  Plot  is  based  on 

illumination  at  <f>  -  w/4. 


forward  direction  lies  in  the  scatterer's  shadow;  that  is,  £lnc  «  -  Escat  in 
the  forward  direction. 

Figure  31  Is  a  four-way  dB  comparison  of  the  monostatic  RCS  computed  by 
both  techniques  for  the  bare  and  coated  magnetic  cylinder.  It  may  be  seen 
that  the  coating  reduces  the  cylinder's  RCS  by  about  13  dB  over  the  fre¬ 
quencies  of  study. 

The  fifth  and  sixth  comparisons  were  the  same  as  the  third  and  fourth, 
but  for  a  perfectly  electrically  conducting  rod.  All  other  parameters  were 
unchanged  between  the  two  pairs  of  comparisons.  Figure  32  gives  the  four¬ 
way  dB  comparison  of  the  monostatic  RCS  for  the  bare  and  coated  electrically 
conducting  cylinder.  It  may  again  be  seen  that  the  coating  results  in  about 
a  13  dB  reduction  of  the  RCS. 
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Scattering  from  a  Thin  Dielectric  Strip 

Richmond^-  has  worked  out  a  solution  for  scattering  of  a  plane  wave  by 
an  infinite  dielectric  strip.  The  strip  may  be  of  any  width,  but  its  opti¬ 
cal  thickness  should  not  exceed  a  tenth  of  a  wavelength. 

This  approximate  solution  is  based  on  superimposing  arbitrary  combina¬ 
tions  of  the  incident  wave  and  the  two  zero-order  waveguide  modes  a  strip 
can  support  propagating  in  its  width  direction.  The  coefficients  of  the 
three  waves  are  then  optimized  by  a  slight  variation  on  Galerkin's  method. 

It  is  an  intrinsic  property  of  the  Richmond  solution  that  the  total 
scattered  current  (conduction,  displacement  and  Prony)  is  symmetric  with 
respect  to  the  center  line  of  the  strip's  thickness  dimension.  Thus,  bis¬ 
tatic  cross-sections  computed  from  this  algorithm  will  always  have  a  plane 
of  symmetry  coincident  with  the  plane  of  the  strip. 

Two  code -code  comparisons  have  been  run  between  the  Richmond 
frequency- domain  algorithm  and  our  TDFD  code.  We  first  did  a  strip  of 
e  -  2  which  was  .04  m  thick  and  2  m  wide.  Cells  .02  m  square  were  used,  so 

1-TiC 

the  TDFD  model  was  2  x  100  cells.  Illumination  at  <f>  -  45 D  with  the 
Poynting  vector  in  the  first  quadrant  was  assumed.  Monostatic  RCS  com¬ 
parisons  were  performed  from  50  MHz  to  500  MHz,  and  bistatic  RCS  comparisons 
were  performed  at  250  MHz  and  500  MHz. 

Figure  33  shows  the  TDFD  model.  Figure  34  Is  a  linear  comparison  of 
the  monostatic  result,  and  Figure  35  is  a  dB  comparison.  Figure  36  is  a 
linear  comparison  of  the  bistatic  result  at  250  MHz,  and  Figure  37  is  a  dB 
comparison.  Figures  38  and  39  are  bistatic  comparisons  at  500  MHz. 

The  second  comparison  was  a  strip  of  -  9  with  ..02  ra  thickness  and 
2  m  width..  Cells  were  again  .02  m  square,  so  the  TDFD  model  was  1  x  100 
c  ills.  All  other  parameters  were  the  same  as  for  the  first  example.. 
Figures  40  -  46  compare  the  same  data  for  this  case  as  Figures  33  -  39  did 
for  the  first  example. 


i 
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illumination  at  <f> 


Dielectric  Strip  Bistatic  Radar  Cross  Section 

0.255  -| 


Bistatic  RCS  for  the  er  -  2  dielectric  strip  at 

250  MHz.  Solid  curve  is  the  Richmond  result, 
and  dashed  line  is  the  TDFD  result.  Plot  is 


Monostatic  Radar  Cross  Section 


Monostatic  RCS  for  the  er  -  9  dielectric  strip. 

Solid  curve  is  the  Richmond  result,  and  dashed 
line  is  the  TDFD  result.  Plot  is  based  on 


illumination  at  £  —  jr/4. 


Bistatic  RCS  for  the  er  -  9  dielectric  strip  at 

250  MHz.  Solid  curve  is  the  Richmond  result, 
and  dashed  line  is  the  TDFD  result.  Plot  is 

1.TIC 

based  on  illumination  at  ^  -  w/4. 


Bistatic  RGS  for  the  «r  “  9  dielectric  strip  at 

500  MHz.  Solid  curve  is  the  Richmond  result, 
and  dashed  line  is  the  TDFD  result.  Plot  is 

3Lt\c 

based  on  illumination  at  <f>  -  w/4. 


It  may  again  be  seen  that  the  two  results  are  in  excellent  agreement 

for  all  frequencies  in  the  —  2  case,  but  that  some  discrepancies  occur 

above  300  MHz  in  the  e  *»  9  case.  For  instance,  the  e  -  9  bistatic  com- 

r  r 

parison  at  500  MHz  (Figure  46)  again  does  not  show  particularly  good 
agreement,  especially  at  the  minor  lobes.  To  get  really  good  finite- 
difference  results,  especially  away  from  field  maxima,  we  again  see  that  a 
wavelength  must  be  resolved  into  16  segments. 

The  TDFD  and  frequency  domain  solutions  are  probably  less  discrepant 
here  in  Figure  46  than  they  were  for  the  er  “  9  dielectric  rod  at  500  MHz 
(Figure  15)  because  we  are  here  only  violating  the  16  cells  per  wavelength 
rule  in  one  direction.  Figure  15  is  based  on  a  calculation  where  the  rule 
w as  violated  in  two  directions . 
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Elliptic  Cylinder  Scattering 


The  transformation  from  Cartesian  to  elliptic  cylinder  coordinates  is 


x  —  1/2  d  cosh  u  cos  v 
y  -  1/2  d  sinh  u  sin  v 
z  -  z 


(1) 

(2) 

(3) 


where  d  is  the  distance  between  the  foci  of  the  ellipse.  Additionally,  we 
sometimes  denote 


£  -  cosh  u 
1 1  —  cos  v 


(4) 

(5) 


so  that 


sinh  u  -  Ji2  -  1  (6) 

2 

In  this  section,  we  shall  use  the  notation  of  Uslenghi  and  Zitran, 

3 

that  of  Blanch  where  it  does  not  conflict  with  the  first  notation,  and  that 
of  Stratton  where  it  does  not  conflict  with  either  of  the  above.  (Lack  of 
a  standardized  notation  greatly  compounds  the  inherently  difficult  problems 
associated  with  elliptic  cylinder  coordinates . ) 


For  all  cylindrical  coordinate  systems,  the  TM  solution  for  Hz  obeys 

V2H  +  k2H  -  0  (7) 

z  z  ' 


In  elliptic  cylinder  coordinates,  this  equation  takes  the  form 


_ 1 

(d/2)2 (cosh2u 


cos2v) 


a2h  a2H  1  a2H 

au2  +  av2J  +  ez2  +  k2Hz 


o 


(8) 


Let  us  now  assume  this  equation  can  be  solved  by  separation  of  vari¬ 
ables  , 
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(9) 


Hz  -  U(u)V(v)Z(z) 
Substitution  of  (9)  in  (8)  yields 


_ 1 _  v«  2^1 

(d/2)2(cosh2u  -  cos2v)  [u  +  Vj+Z  +  ^  "*  ^ 

Let  C  be  the  first  separation  constant, 

7." 

z  +  c  -  0 


Then  U  and  V  must  obey 


U"  -  (a  -  1/8  •  d2(k2  -  C)  cosh  2u)U  -  0 


V"  +  (a  +  1/8  •  d2(k2  -  C)  cos  2v)V  -  0 


(10) 


(ID 


(12) 


(13) 


where  a  is  the  second  separation  constant. 

In  the  rest  of  this  report,  we  shall  assume  there  is  no  z  dependence, 
so 


k2  -  C  -*  k2  (14) 


and  we  shall  denote 


c  -  1/2  •  (kd)  (15) 
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where  c  is  not  the  speed  of  light,  but  is  the  number  of  radians  at  the 
frequency  of  interest  between  the  ellipse  center  and  one  of  its  foci.  Then 
eqs.  (12)  and  (13)  become 

U"  -  (a  -  1/2  •  c2  cosh  2u)U  -  0  (16) 

V"  +  (a  +  1/2  •  c2  cos  2u)V  -  0  (17) 

Equation  (17)  is  Mathieu's  equation  and  (16)  is  the  modified  Mathieu 
equation. 

The  c  of  eqs.  (15)  -  (17)  corresponds  to  c  \  in  Stratton  and  Js  in 

°  5 

Blanch.  In  some  more  modern  works,  such  as  Abramowitz  and  Stegun,  and 
Hodge, ^  it  corresponds  to  2/q. 

It  turns  out  that  Mathieu's  equation  only  permits  periodic  solutions 
for  discrete  eigenvalues  of  a,  where  these  eigenvalues  depend  on  c.  If  we 
assume 


V 


2  De9^(c)  cos  2kv 

k-0 


(18) 


we  obtain  even  solutions  of  period  tt.  The  rth  eigenvalue  of  this  equation 
2r 

is  denoted  a  ,  and  the  rth  eigenfunction  is  denoted 


V 


CO 


I 


Be 


2r 

2k 


(c)  cos  2kv 


(19) 
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If  we  normalize  these  functions  to  be  unity  at  v  - 
Mathieu  functions  of  period  ir 


V  2r 

Se2r(c,»7)  -  )  De2k(c)  cos  2kv 

k=0 


Alternatively,  if  we  let 


V  “  )  De2k+l(c)  cos<21:  +  1>v 
k-0 


we  obtain  the  even  solutions  of  period  2jt.  The  rth 
system,  also  normalized  to  be  uni.ty  at  v  -  0,  is 


Se2r+l^C’'7'  ~  2  De2k+l(c^  cos<2k  +  1)v 
k-0 


Now  let  us  consider  odd  solutions  of  period  jt, 


V 


2  Do2k(s)  sin  2k- 
k~l 


2r 

The  rth  eigenvalue  of  this,  equation  is  denoted  b 
tion,  normalized  no  have  unity  derivative  at  v  «  0  is 


0,  we  obtain  the  even 


(20) 


(21) 


eigenfunction  of  this 


(22) 


(23) 


and  the  rth  eigenfunc- 
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SoZr(c,»7) 


-  5 


Doz^(e)  sin  2kv 


k-1 


Finally,  the  rth  odd  eigenfunction  of  period  2tt  is 


(24) 


*o 

5*, 

S5 


So2r+l(c,,?J  "  2  D°2k+1^  sin^2k  +  1>V 
k-0 


(23) 


3i^ 

I 


Any  even  Mathieu  function  is  orthogonal  to  any  odd  Mathieu  function  on 
the  interval  (0,2ir).  Likewise,  different  even  or  odd  Mathieu  functions  are 
orthogonal,  The  normalization  factors  are 


2r 

N^(c)  “  ^  Se2r(c,cos  v)2  dv 


?r 


2(De*r(c))2 


(De*r(c))2 


-] 


(26) 


(c> 

w2r+rc' 


*  )  {De2k+l(c))S 
k-0 


(27) 


N^(c) 


*  Y  0>°2>!<c)>2 


(28) 


^;i<c> 


V  2r+l ,  . . 

*  /  trVwl(c)) 

k-»0 


(29) 
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The  four  types  of  Mathieu  functions  defined  by  eqs.  (17)  -  (29) 
describe  the  azimuthal  dependence  of  the  elliptic-cylinder  harmonics. 
Equation  (16),  on  the  other  hand,  describes  the  radial  dependence.  There 
are  a  total  of  16  kinds  of  radial  Mathieu  functions. 

Let  Z^^(x)  denote  the  kth  Bessel  function  of  the  jth  kind.  For  ex¬ 
ample  , 


z£3)(x>  -  Jk(x)  +  iYk(x) 


(30) 
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Then  the  four  radial  Mathieu  functions  corresponding  to  Se2r(c,»j)  may  be 
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shown  ’  to  be 
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where 
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and  s  is  any  arbitrary  interior,  but  is  best  selected  if 
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Similarly,  the  four  radial  Mathieu  functions  corresponding  to 
Se2 r+i_(c,9)  are 
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The  radial  Matnieu  functions  corresponding  to  So2r+p(c,»j)  are 
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Let  us  now  consider  the  expansion  of  a  plane  wave  in  elliptic-cylinder 

inc 

harmonics.  If  we  have  a  plane  wave  propagating  at  an  angle  <f>  from  the 
+  x  axis, 
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then  the  desired  expansion  is 
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The  associated  electric  field  is 


,inc  .  - 1 
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where  h  is  the  elliptic  cylinder  metric 


h  -  hy  -  hv  -  (d/2)(cosh2u  -  cos2v)J 
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Expansion  of  eq.  (43)  yields 
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We  shall  first  assume  this  plane  wave  is  scattered  by  a  perfectly 
conducting  elliptic  cylinder.  Thus,  the  scattered  solution  must  be  an 
infinite  series  of  outgoing  Mathieu-Hankel  functions. 
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The  requirement  that  on  {  «  (j 
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We  are  now  in  possession  of  enough  information  to  evaluate  the  RCS  of 

(3)  (3) 

the  elliptic  cylinder.  At  large  arguments,  Re  (c,£)  and  Ro  (c,£)  both 
,4  mm 

approach 
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Uslenghi  an  Zitran  then  define  the  scattering  function 
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Vjrkr* 


i(kr'  -  «/4) 
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where  approaches  <p ^  at  large  In  view  of  eqs.  (I),  (2),  and  (15),  we 

can  show  that 


kr'  -♦  c  cosh  u 
P 


(53) 


Substitution  of  (46),  (49)  -  (51),  and  (53)  into  (52)  then  yields  the  bis- 
tatic  scattering  function. 
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We  can  then  evaluate  the  bistatic  RCS, 


RCS(Vp) 
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Uslenghi  and  Zitran  have  published  curves  for  |P(v^)|2  as  functions  of 
,  and  tanh  ux .  The  quantity 


c$i  -  (1/2)  •  kd  cosh  Ui 


(56) 


is  the  number  of  free  space  radians  at  the  frequency  of  interest  along  one 
half  of  the  major  axis  of  the  ellipse. 

Additionally , 


tanh  ux  -/$!  -  l/$i 


(57) 


is  zero  if  the  ellipse  is  reduced  to  a  conducting  strip,  and  is  unity  in  the 
special  case  where  the  ellipse  fattens  to  a  circular  cylinder. 


Figures  47  -  49  illustrate  our  results  for  |P(v')|2  when  $inc  is  0,jt/4 

and  w/2  and  v'  is  <f>  +  x  (the  backscatter  case).  These  figures  agree 

P  2 

exactly  with  previously  published  results  and  confirm  tne  correctness  of 
our  analysis  and  Mathieu  function  routines. 


We  have  also  reproduced  the  Uslenghi  and  Zitran  results  for  TE  scat¬ 
rering  by  an  elliptic  cylinder,  although  this  calculation  is  only  relevant 
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to  our  present  work  in  that  it  tests  our  Mathieu  function  routines  under 
different  circumstances. 

It  is  interesting  to  note  that  in  the  case  of  scatterii  g  by  a  conduct- 

(IV  scat 

ing  strip,  ux  is  zero,  £x  is  unity,  Re^  (c,l)  is  zero,  and  all  the  am 
vanish. 

For  scattering  by  a  dielectric  cylinder,  we  encounter  a  strange  com¬ 
plication  which  will  probably  be  new  to  anyone  who  has  previously  only  used 
Bessel  functions  and  spherical  harmonics.  The  incident  and  scattered  fields 
can  still  be  represented  by  eqs.  (40)  -  (47).  However,  we  must  now  also 
consider  the  fields  which  penetrate  the  cylinder, 
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Here  YA  is  the  characteristic  admittance  of  the  dielectric 


Yi  -/e i//*c 


kx  is  the  wavenumber  in  the  dielectric 


ki  “  “V^o 


and  Cj  is  the  number  of  radians  between  the  ellipse  center  and  a  foci  for  a 
wave  travelling  in  the  dielectric. 


cx  -  (1/2)  •  ktd 


The  boundary  conditions  at  f  -  |x  are  now  the  continuity  of  H  and  E  , 
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These  conditions  lead  to  the  relationships 
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The  complication  now  arising  is  that  Se  (0,17)  and  Se  (c ,,>j)  are  not 

m  u\ 

equal,  since  c  and  cx  are  different.  ';..u  other  words,  the  dielectric 
cylinder  c  rnses  each  incident  elliptic  harmonic  to  couple  to  every  scattered 
elliptic  harmonic.  Thus,  unlike  the  case  for  a  circular  cylinder  or  a 
sphere,  eqs  (65)  and  (66)  do  not  separate  out  into  decoupled  equation  pairs 
fcr  each  value  of  m. 

It  is  necessary  to  substitute  eqs.  (20)  and  (22)  for  Se  (c,q)  and 
Sem(cj ,ij)  into  (65)  and  (66).  Doing  thi3,  rearranging  terms,  and  factoring 
out  the  trigonometric  functions  yields 
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These  equations  must  be  satisfied  for  all  n,  as  must  a  similar  family  of 
equations  for  the  odd  elliptic  harmonic-  Once  such  a  solution  has  been 
performed,  the  RCS  is  still  given  by  eqs.  (AS)  -  (53).  We  have  not  yet 
programmed  this  solution,  but  certainly  have  all  the  necessary  bits  and 
pieces  of  Mathieu  programs  on  the  shelf,  and  could  quickly  assemble  them. 
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The  two-medium  elliptic  cylinder  problem  is  a  fairly  straightforward 
extension  of  eqs.  (67)  and  (68).  Unlike  Eessel  functions  of  the  second 
kind,  radial  Mathieu  functions  of  the  second  kind  do  not  go  to  infinity  at 
u  -  0.  However,  they  do  have  both  nonzero  values  and  nonzero  derivatives 
there.  This  means  elliptic -cylinder  harmonics  with  radial  Mathieu  functions 
of  the  second  kind  have  either  a  sharp  ri<ge  or  a  step  discontinuity  on  the 
line  connecting  the  coordinate  system  foci.  Consequently,  their  appearance 
is  forbidden  in  the  case  of  scattering  from  a  uniform  single-medium  elliptic 
cylinder.  Such  is  not  the  case,  however,  for  the  outer  medium  in  the  two- 
medium  case. 

Subject  to  this  understanding,  the  elliptical -cylinder  harmonic  coeffi¬ 
cients  for  the  two-medium  cylinder  obey 
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The  odd  elliptic-cylinder  harmonic  coefficients  obey  a  similar  equation 
family . 

If  the  inner  region  is  a  perfect  conductor,  we  need  only  be  concerned 

with  the  vanishing  of  the  total  tangential  electric  field  at  its  surface. 

Thus,  in  that  case,  we  can  ignore  eq,  (69)  which  pertains  to  continuity  of 

H  .  rforc'wer,  we  know  that  a^nner  will  vanish  if  the  inner  material  is  a 
z  m 

perfect  conductor.  This  has  ul.e  effect  of  decoupling  the  elliptic  harmonics 
from  each  other  at  the  inner  interface;  eq.  (70)  reduces  to 


rvi^w  +  *rer<2)  ^i2>'(=,.«.>  -  o 


Equations  (71)  and  (72)  do  not  simplify. 

The  problem  of  a  perfectly  conducting  ellipse  or  strip  confocaily 
coated  by  a  dielectric  is  probably  by  far  the  closest  replication  of  a 
stealth  wing  wbicYi  there  is  any  chance  of  solving  canonically.  As  such, 
this  problem  is  cf  gre^c  importance  in  testing  out  numerical  codes  for  RCS 
reduction. 
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SUMMARY 

In  this  report,  we  have  presented  canonical  solutions  for  2D  TM 
electromagnetic  scattering  by  lossy,  layered  circular  cylinders;  thin 
dielectric  strips;  thin  perfectly  conducting  strips;  and  lossless,  confo- 
cally  layered  elliptic  cylinders. 

These  canonical  solutions  have  been  used  to  check  out  a  TDFD  code 
designed  to  give  the  monostatic  and  bistatic  RCS  of  a  generalized  cylinder. 
(The  TDFD  result  undergoes  a  Fourier  and  a  near-field  to  far-field  transfor¬ 
mation  before  yielding  an  RCS  which  can  be  checked  against  the  canonical 
results . ) 

The  TDTD  code  can  handle  cylinders  containing  abrupt  electrical  discon¬ 
tinuities,  including  conducting  or  resistive  cards,  anisotropy  in  the  plane 
of  the  cylinder,  and  frequency  dependence  in  e,  fi  and  a. 

The  only  major  limitation  we  have  found  on  the  TDFD  code  is  tnat  about 
16  cells  are  required  to  resolve  the  shortest  wavelength  of  interest.  Also, 
about  80  cells  must  be  interposed  between  the  scatterer  and  the  outer  bound¬ 
ary  to  isolate  the  reactive  field  of  the  scatterer  from  outer  boundary 
effects. 
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APPENDIX  1 

TIME-DOMAIN  TREATMENT  OF  MAXWELL'S  EQUATIONS 
IN  FREQUENCY  DEPENDENT  MEDIA 


Introduction 

Consider  a  medium  with  anisotropic,  frequency- dependent  electrical 
properties.  The  electrical  response  of  such  a  material  may  be  fairly 
generally  described  by 


V  x  H  -  J(t)  +  Jf(t) 


where  J^Ct)  is  a  forced  current  and 


1(b)  “  £o  *  E(t)  +  ’  DE(t)  +  K(t-t')  *  DE(t')dt’ 


with  D  indicating  the  time-derivative  operator.  In  this  formulation,  o0, 
e  ,  and  K(t-t')  are  second-rank  tensors. 

mO  '  im  '  " 


The  frequency- domain  form  of  eq.  (2)  is 


J(u)  «■  o0  +  iwe  + 


iw  J^e  iwUK(u)duj  *  E(w) 


Separation  of  eq.  (3)  into  real  and  imaginary  parts  gives  representations 
for  the  frequency- dependent  conductivity  and  permittivity  tensors, 


£(«)  -  2,0  +  Ee  j^iwj^e "  ^WUK (u)  duj 
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(9) 


Longmire  and  Longley  assuaed  that  materials  could  be  represented  by  the 
exponential  series  of  eq.  (6)  yith  one  term  for  each  decade  of  frequency 
over  the  spectrum  of  interest,  'inis  is  equivalent  to  doing  a.  Prony  expan¬ 
sion  of  K(u)  [or  a(u)  and  e(w)]  with  the  poles  forced  to  be  spaced  at 
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While  this  assumption  has  been  claimed  to  be  reasonably  accurate  for  wet 
soil,  it  would  seem  generally  more  correct  to  determine  the  poles  from  a 
Prony  analysis  of  the  medium's  measured  frequency- dependent  characteristics. 
This  is  especially  likely  to  be  true  if  the  material  exhibits  rapid  varia¬ 
tion  in  a  and  e  with  frequency. 

State  Theory  Applications 

let  us  first  assume  the  Prony  analysis  reveals  no  complex-conjugate 

pole  pairs.  In  general,  the  a  will  be  second  rank  tensors,  but  the  P  will 

**m  'm 

only  be  scalars.  Then  for  every  pole,  each  component  of  J  will  obey 


DE1  •  <Wml  -  Val  -  0  »  -  1  -  M.  t  -  1  -  3 


(11) 


Additionally,  the  tensor  form  cf  eq.  (7)  giv/e- 


H 


€«ijDEj  +  a°ijEj  +  £1amij  Jmj  “  Jj  i,j  «  1  -  3 


+  aoi  .jEj  +  )  a  ..  J  .  -  J 


(12) 


where 


J.  -  (V  x  H  -  Jf)j 


(13) 


Equations  (11)  and  (12)  constitute  a  set  of  3(M+1)  coupled  first  order 
differential  equations . 
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If  anisotropy  and  frequency  dependence  were  not  present,  the  usual 

method  of  numerical  solution  would  be  explicit  time-domain  finite 

differencing.  In  this  method,  E  and  H  evaluation  points  alternate  bf>th 

8  - 10 

spatially  and  temporally  using  a  well-tested  leapfrog  arrangement.  In 

this  arrangement,  no  two  equations  are  coupled,  and  E^+^^(I,J,K)  means  Ex 
evaluated  at  ((I  +  1/2)AX,  JAY,  KAZ,  (n  +  l/2)At) . 

However,  the  present  system  of  equations  requires  the  three  E^'s  and  3M 
Jmj  's  all  to  be  evaluated  simultaneously.  While  this  cannot  be  done  using 
conventional  time-domain  finite  differencing,  state  theory  does  indicate  an 
appropriate  generalization  of  time-domain  finite  differencing. 


First,  let  us  consider  the  case  where  anisotropy,  but  not  frequency 
dependence,  is  present, 


i  s 
$  & 


k 


[ej  D[E]  +  [er0][E]  -  [J] 


This  matrix  differential  equation  has  a  homogeneous  solution 


(14) 
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[Ek  -  e 


-UJ 


[A] 


(15) 
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and  a  particular  solution 
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eJ  Jv1. 
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Nr 


[E]p  -  [cxo]*1!!] 


giving  a  general  solution 


-Uj'Volt  .! 

[E]  -  e  *  [A]  +  [a0]  A[J] 


(16) 


(17) 
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The  constant  vector  [A]  may  be  evaluated  at  (n  -  1/2) At: 


[E]n‘1/2  -  [A]  +  [cr0] 


(18) 


This  gives  the  new  E-field  vector  in  terms  of  the  old, 


[E] 


"1/2  „  e‘l£“1  ll<7°li,:lE!n-l/2  +  ( 


[I]  -  e 


_  i.ltt in 


[»0]  [J1  (19) 


Similar  exponential  matrix  techniques  have  been  reported  for  time- 
domain  solution  of  generalized  multi-conductor  transmission  lines. ^  In  the 
previous  work,  one  may  see  how  to  evaluate  eq.  (19)  if  [ct0]  is  singular  or 
if  [e  ] ~^[cr0]At  has  arbitrarily  large  elements.  Basically,  matrices  are 
exponentiated  using  the  power-series  representation  of  an  exponential. 

If  frequency  dependence  is  present,  the  [E]  vector  of  eqs.  (14)  -  (19) 
becomes  replaced  by 


[E] 


(20) 


The  f«  1  matrix  becomes 

1  «oJ 


and  the  [o0]  matrix  becomes 


(21) 
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Lastly,  the  forcing  vector  becomes 


[J]  ->  0  -  [J'] 


Then  the  matrix  equation  for  simultaneous  advancement  of  E  and  the  is 


^n+1/2  _  e-[«'l  [o' lAt[E» ]n-l/2  +  [[j]  .  e" t£f 1  ia' ]Atj [a, ]n 


*  1  /»» 


In  the  past,  time-domain  finite  differencing  has  not  often  considered 
anisotropy.  Frequency -dependent  effects  have  been  included  by  using  the  old 
to  find  the  new  £n+^^.  (This  decouples  £  from  the  In  eq.  (12).) 
Then  the  new  have  been  used  to  find  the  new  fr0m  eq.  (11). 


Treatment  of  Complex  Poles 

If  Prony  analysis  of  the  material's  frequency  dependence  reveals  com¬ 
plex  pole  pairs,  a  more  general  treatment  becomes  necessary.  In  this  case, 
K(u)  will  contain  terms  of  the  form 
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b  sin(7  u  +  <f>  )e 
«m  m  m 


-  b  co s<i  sin7  ue  +  b  sin A  COS7  ue 
-m  rm  m  -m  m  'm 


The  J  (t)  of  eq.  (8)  now  becomes 


J^/t)  -  J  DE(t')(cos^msin7m(t-t')e  m 


-j9 

+  sin^mcos7m(t-t' )e  m  )dt' 


-  .J  ( t)cosi  +  J  (t)sini 
“me  m  Tns  m 


where  J  and  J  are  the  parts  of  J  associated  with  cosA  and  sin^  , 
“me  “ms  “m  m  m 

respectively: 


^(t)  -  J  DE(t')sin7m(t-t')e  dt 

-00 


,(t')  -  J  DE(t')cos7m(t-t') 


e  dt' 


Differentiation  of  J  (t)  and  J  (t)  yields 

“me  tus 


DJ  (t)  -  -5  J  (t)  +  7  J  (t) 
^acv  1  P7sttslck  '  Tims'-  J 
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(t)  -  DE(t) 


0  J  (t) 
rm“nis '  ' 


+  7  J  (t) 
'nnmc'  7 


(30) 


These  equations  can  be  solved  for  J  (t)  and  J  (t)  using  Keavyside  algebra: 


i<D  +  W  +  4c<t>  -  yj*™  <31> 

i<D  +  Vs  *  ■>%  ^(w  -  <“  +  ejme>  .  <32> 


Thus,  J  (t)  obeys  the  differential  equation 


[CD  +  euV  +  7j]  V t)  -  [iln#.CD  +  +  coS(Sm7m]DI(t)  (33) 


In  principle,  equations  like  (33)  could  be  added  to  the  set  of  equa¬ 
tions  given  by  (11)  and  (12),  and  the  entire  ensemble  solved  by  state 
theory.  This  approach,  however,  requires  treatment  of  second- order  matrix 
differential  equations  of  the  form 


[A]D*[E]  +  [B]D(E]  +  [Cl [E]  -  [F] 


(34) 
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The  homogeneous  solution  of  this  equations  includes  square  roots  and  complex 
exponents  of  matrices;  it  is  much  more  complicated  than  eqs.  (14)  -  (19). 
(To  the  best  of  our  knowledge,  exponential  differencing  has  never  been 
applied  even  to  scalar  second-order  differential  equations.) 

Consequently,  when  Prony  analysis  of  the  material  data  yields  complex 
pole  pairs,  our  present  strategy  is  to  fall  back  to  the  old  technique  for 
dealing  with  real  poles:  First  find  En+V2  using  the  old  ,2^  .  Then  use 
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the  new  En+1</2 
Tn+l/2 

T1 


end  the  finite-difference  form  of  eq. 


(33)  to  find  tha  new 


» 


i 


121 


